Water slippage on non-ideal surfaces: the role of static and dynamic roughness 
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We present non-equilibrium Molecular Dynamics simulations of water flow in nano-channels with 
the aim of discriminating static from dynamic contributions of the wall atoms to the slip length 
of the molecular flow. To this purpose, we vary the rigidity of model springs and the mass of the 
tethered wall molecules independently, so as to separately tune the rugosity and the flexibility of 
the walls. First, it is found that the slip length for small to moderate values of the shear intensity is 
significantly enhanced by increasing the wall flexibility. Second, and perhaps more importantly, we 
show that the regularization of the slip length divergence at high shear rates, formerly attributed 
to wall flexibility, is controlled only by the static surface roughness. Minimal displacements, caused 
for example by thermal motion at room temperature, are shown to be sufficient to remove the 
aforementioned divergence. 

PACS numbers: 47.11. Mn,47.61.-k,83.50.Rp 
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In the pioneering paper [30(, it was shown that the 
Navier slip boundary condition u s = (. s ^ can be re- 
garded as the low-shear limit of a universal, nonlinear 
relation between slip velocity u s , the slip length £ s , and 
the local shear rate 7. In the approximately fifteen years 
since then, violations 3ll433||. as well as confirmations 
34M38J] of this behavior have been reported, preventing 
the emergence of a clear consensus on the dependence of 
fluid slippage properties on cither the imposed shear, or 
the static and dynamic surface roughness. This lack of 
consensus calls for further investigations, since control- 
ling and predicting the behavior of water in nanochan- 
nels, where surface properties have a profound effect on 
flow dynamics, stands as a major scientific challenge with 
many practical applications in modern science and tech- 
nology. Indeed, an accurate understanding of nanoscale 
friction phenomena at fluid-solid interfaces is paramount 
to the design of micro and nanofluidic devices aimed at 
optimizing mass transport against an overwhelming dis- 
sipation barrier [39h41| . The effects of static properties 
of the fluid-solid interface, such as the contact angle and 
surface roughness, on slippage phenomena have been in- 
vestigated in depth in the recent past |30l. I33J 1371. \38L 142- 
l46j . However, the role of dynamic properties, such as the 
wall flexibility and thermal conductivity, remains much 
less explored [32|,|47( . Most importantly, the geometric ef- 
fects due to corrugations of flexible walls, have never been 
disentangled from the dynamic ones. A recent analysis 
pointed out that all investigations to date which report 
absence of a divergent slip-length under large shear were 
performed with flexible walls capable of exchanging mo- 
mentum and energy with fluid [481 ] . From this observa- 
tion, and supported by additional numerical calculations 
and analytical models, the authors of Ref. [48j concluded 
that heat and momentum transfer from the fluid to the 
wall are responsible for the observed saturation of the 
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FIG. 1. Simulation snapshots of a portion of the solid sur- 
face and corresponding potential energy isosurfaces for a sin- 
gle water molecule at energy ksT , 2ksT and 3ksT, respec- 
tively from the highest to the lowest color saturation. The 
vertical position of the isosurfaces is magnified by a factor 
5. Panel (a): standard surface with no atomic displacement; 
Panel (b) : random quenched functionalization with no atomic 
displacement; Panel (c): random quenched functionalization 
with atomic displacement £ = y/kbT '/k ~ 0.022 nm. 



slip length at high shear rates, as opposed to the diver- 
gent behavior which takes place using rigid walls. This 
picture has been questioned in [49(, where the authors 
attributed the divergence to the use of a thermostat ap- 
plied to the fluid molecules. Instead of a saturation of 
the slip length at high shear rate, these authors report a 
slippage drop toward the non-slip condition due to fluid 
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FIG. 2. Slip length as a function of the shear rate for different 
surface types. Very small values of excursion £ remove the slip 
length divergence for both surface types. Wall atom mass 
m = 1.2 x 10 3 amu. 



heating. With the amount of slippage at high shear rates 
being reported by different authors to cither vanish, reach 
a constant finite value, or to diverge, the present picture 
appears rather controversial. More importantly, no clear 
consensus has yet emerged on the role of different physi- 
cal quantities governing the transition from a finite to a 
divergent slip length. 

We present here the results of a series of simulations 
performed with the aim of discriminating those contri- 
butions to the slip length which are of purely dynamical 
origin, from those whose origin is rooted in the surface 
configurational properties. To this end, we have per- 
formed simulations of water flow in channels with flexible 
walls varying the masses to of the wall atoms and, inde- 
pendently, the harmonic constant k of the potential that 
bounds wall atoms to their lattice sites (see Fig. [IJ. Ac- 
cording to our findings, this approach proves key to shed 
new light on the slip length behavior at high shear rates. 
More precisely, we report numerical evidence that it is 
the presence of an even extremely tiny amount of disor- 
der in the atomic positions at the surface, rather than 
wall flexibility, which proves responsible for taming slip 
length divergence and making it shear-independent. 

Using the Gromacs suite [50|, we simulated a slab 
of water molecules confined between two graphitc-likc 
walls, each consisting of three atomic layers. The water 
molecules are modelled using the SPC/E potential [511 ] . 
and the electrostatic interaction is calculated using the 
smooth particle mesh Ewald method [52[, together with 
Yeh-Berkowitz correction 53j to remove the contribution 
from periodic images along the wall normal. The distance 
between the two atomic planes in direct contact with wa- 
ter defines the channel width L = 6.306 nm, while the 
entire, rectangular simulation cell is 13.035 by 16.188 by 
8.0 nm along x, y and z, respectively, with the wall sur- 
face normals along the z-direction. The integration time 



step for the equations of motion was 1 fs. The slip length 
i s is calculated based on its its geometrical definition 



t a = n/X - L/2, 



(1) 



where the viscosity n 
f/u v 



A 



41 



f/j and the friction coefficient 
are determined by the force per unit 
surface between wall and water, /, the wall speed u w 
and the shear rate 7. 

We investigated two different setups, with respect to 
the wall-water interatomic interaction potential. In the 
first setup, hereafter identified as "standard", the inter- 
action between wall and the oxygen atoms in the wa- 
ter molecules is prescribed by a Lennard- Jones poten- 
tial U(r) = C 12 r- 12 - C 6 r- 6 with C 6 = 2.47512 x 10~ 3 
kJ/nm 6 and C 12 = 2.836328 x 10" 6 kJ/nm 12 [13 up to 
an interatomic distance r of 0.9 nm. Above that distance, 
the potential is smoothly switched to zero at r — 1.2 nm, 
using a fourth order polynomial. The second setup is a 
random quenched functionalization of the first, realized 
by making 40% of wall atoms purely repulsive (Cg = 0), 
and increasing the interaction strength of the remaining 
ones by a factor a = 1.977. This generates an intrinsic, 
mainly horizontal rugosity of the surface: see panel (b) 
in Fig.[TJ The value of a has been chosen such that water 
equally wets the surfaces of the two setups, resulting in 
comparable macroscopic contact angles [55|. Given the 
known dependence of slip length on contact angle [45(, 
equal wetting is prerequisite to comparing the slip length 
along two microscopically different surfaces and under- 
standing to what extent the slip length is affected by dif- 
ferent types of surface inhomogeneities, rather than by 
the contact angle. 

For each of the surface type just described, we simu- 
lated a flexible and a rigid variant. The flexible variant 
is realized by tethering the atoms of each wall's outer- 
most layer to their respective lattice sites via a harmonic 
potential Uk(r) — \kr 2 . The spring constant k defines 
the characteristic excursion £ = ^JksT/k of the teth- 
ered atoms, so that a small value of k leads to a high 
effective rugosity. The rigid variant, on the contrary, is 
realized by freezing all surface atoms either at their lat- 
tice site position (corresponding effectively to £ = 0), or 
in a configuration taken from the equilibrium trajectory 
of the flexible variant with finite £, which wc will denote 
by to = 00. Both mobile and fixed wall atoms interact 
with water only, and not among themselves. By varying 
the mass to of the tethered molecules, it is possible to 
achieve a separate control over dynamic properties, as 
the value of to does not influence configurational prop- 
erties (e.g.: particles distribution, free energies, contact 
angles), but only dynamical ones, like water-surface fric- 
tion and thermal conductivity. As a result, we are able 
to assess the dependence of the slip length on dynamic, 
static and non-equilibrium properties, controlled by the 
mass and mean excursion of the wall atoms and by the 
imposed shear. 
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F IG. 3. Slip length as a function of the excursion £ = 
\JksT ' jk for two different values of shear rate (squares: 7 = 
0.0125 ps~ : ; circles: 7 = 0.625 ps _1 ). Tethered atoms have a 
mass of 1.2 x 10 3 amu. 



Fig. Q] shows the potential energy isosurfaces for three 
different cases. While the spatial functionalization (b) 
already introduces a noticeable perturbation of the en- 
ergy landscape, the addition of an even tiny spatial dis- 
placement of the surface atoms (c) generates the largest 
perturbation. 

As a first step, with the mass of surface atoms held at 
a constant to = 1.2 x 10 3 amu (to efficiently integrate the 
equations of motion of the tethered atoms even with the 
largest spring constant) we investigated the shear rate 
dependence of the slip length for different displacements 
£. Providing useful information about the scales which 
come into play in the determination of slippage proper- 
ties, we present an overview of this dependence for both 
the standard and randomized surfaces in Fig. f2] 

Five different cases were explored for the random 
quenched surface, between the values £ = and £ ~ 0.11 
nm (A; = 200 kJ/mol nm -2 ), the former mimicking an 
infinite spring constant, and the latter corresponding to 
the one proposed as optimal in [48[. The smallest but 
still finite value, £ ~ 0.0027 nm, was chosen to model a 
nearly flat wall which yet retains dynamical features. 

When the wall atoms are perfectly flat (£ = 0), both 
surface types show a tendency towards a larger — in fact 
divergent — slip length at increasing shear. It takes only 
sub-nanoscopic excursions of wall atoms to remove the 
divergence and make the slip length shear independent. 
The slip length is shear-independent (constant) also for 
large values of £, however the constant value does increase 
with the spring constant k (decreasing excursion £) to a 
limiting value of about 8 nm for the random quenched 
surface and a slightly larger one for the standard surface. 
The maximum limiting value for shear-independent slip 
lengths is reached at £ c ~ 0.01 nm. For smaller mean ex- 
cursions, the slip length is no longer shear-independent 
and a divergent-like behavior is observed instead. At first 
glance, a characteristic excursion £ = 0.01 nm might ap- 
pear to be surprisingly small, but we note that several 
other important scales in this system, such as the width 
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FIG. 4. Left axis: The density distribution of tethered atoms 
(thin lines) and water molecules (thick lines) in the random 
quenched case at £ ~ 0.07 nm (solid) and £ ~ 0.016 nm 
(dashed). Right axis: water velocity profiles at shear rate 
7 ~ 0.016PS" 1 (squares: £ ~ 0.07 nm; circles: £ ~ 0.016 nm). 



of the lksT energy basin of a water molecule in the direc- 
tion normal to the surface (see Fig. []}, and the distance 
5 ~ 0.03 nm travelled by a water molecule during its 
translational relaxation time of r r ~ 0.1 ps [561 ]. 

The transition from shear-independent to divergent 
slip length is made even more evident by plotting the 
slip length as a function of £ (Fig. [3]) for two selected 
shear rates, 7 ~ 0.016 ps -1 at the upper end of the low 
shear-rates plateau, and 7 ~ 0.08 ps -1 well into the shear 
region where divergence is observed. At the lower shear 
rate, saturation is reached for excursions smaller than 
£ c ~ 0.01 nm, whereas the high-shear rate slip length 
keeps increasing, seemingly unaffected. The wall and wa- 
ter density profiles are more sharply peaked in the small 
excursion case, but the effect is significant only next to 
the surface, and diminishes greatly by the second density 
peak. In both cases the flow velocity is well represented 
by the Couette flow solution, even in regions near the 
wall where the density of water becomes negligible. 

From the analysis presented thus far, the role of wall 
flexibility remains unclear, as this dynamical aspect was 
not separated out from the other features. For this rea- 
son, we selected two extreme values of excursion (£ = 
0.0027 and 0.11 nm) and calculated the slip length for 
several values of the tethered masses to, including the 
case of atoms fixed in one equilibrium configuration, 
to = 00, as described before. 

Results from these calculations are shown in Fig.[SJ At 
low shear rates (7 < 10~ 2 ) both values of £ demonstrate 
marked changes to the slip length with the mass m. The 
small excursion case £ = 0.0027 nm halves in slip from 
the lowest mass to = 1.2 x 10 3 amu to to = 00, and in 
the £ = 0.11 nm case the slip length is reduced to zero at 
to = 00. The transition from thermally insulating walls 
(to = 00) to conducting walls is therefore associated to 
a substantial increase in slip length. In other words a 
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FIG. 5. Slip length as a function of the shear rate in the 
random quenched case, for £ ~ 0.0027 nm (squares: m — 
1.2 x 10 3 uma; circles: m = 1.2 x 10 5 uma; upper triangles: 
m = oo) and for £ ~ 0.11 nm (lower triangles: m = 1.2 x 
10 3 uma; rhombs: m = oo). Lower panel: water viscosity as a 
function of shear rate. The experimental value ffcxp = 515.5 
kj/mol/ps/nm 3 at 300K has been interpolated from reference 
data[57(|. 



flexible wall reduces friction on the fluid by being capable 
of adjusting locally to the flux of water molecules. In the 
extreme case of m = oo, the slip- length divergence is 
indeed removed, provided that the excursion £ is larger 
than £ c . This provides evidence that the flexibility can 
not be responsible for the regularization of the slip length 
at high shear-rates. 

It is important to note that these are entirely surface- 
related effects. The viscosity of water r\ is the only bulk 
property that appears in the definition of the slip length 
for a Couette flow, Eq. ((T|), and it does not show any 
dependence on shear, or spring constant, as shown in 
Fig. [5j bottom panel. Therefore, we argue that the tran- 
sition from a divergent to a shear-independent slip length 
is determined not by the wall flexibility, but rather by the 
configurational disorder induced by tiny displacements of 
the wall atoms. To get an appreciation for the physical 
realism of the characteristic excursion lengths discussed 
here, we compare them to displacements obtained from 
the effective spring constants of real graphite at 300 K 
measured by inelastic X-ray scattering [58[ , which corre- 
spond to about 0.006 nm for nearest neighbors, and are 
much larger for second neighbors. Thermal fluctuations 
alone are therefore expected to be sufficient to attain a 
shear-independent slip length. 

Summarizing, we have performed non-equilibrium 
Molecular Dynamics simulations of water flow in nano- 
channels with separate control of the flexibility and 
roughness of the wall. The simulations show that the 
disappearance of the divergence of the slip length at high 
shear rates, formerly ascribed to wall flexibility, is due 
instead to the excursion of wall atoms from the ideal hor- 



izontal plane. Moreover, these results show that atomic 
displacements as small as those due to thermal motion 
at room temperature are sufficient to regulate the diver- 
gence of the slip length. We conclude from this that, even 
in the absence of surface imperfections such as point de- 
fects or dislocations, high-shear water flows in nanoscale 
channels should not exhibit any divergent slip length; 
that thermal motion of the wall atoms is sufficient to 
tame such divergence. 
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